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Fhe dynamics of the three coupled bosonic wells (trimer) containing N bosons is investigated 
within a standard (mean-field) semiclassical picture based on the coherent-state method. Various 
periodic solutions (configured as 7r-like, dimerlike and vortex states) representing collective modes 
are obtained analitically when the fixed points of trimer dynamics are identified on the iV=const 
submanifold in the phase space. Hyperbolic, maximum and minimum points are recognized in the 
fixed-point set by studying the Hessian signature of the trimer Hamiltonian. 

Fhe system dynamics in the neighbourhood of periodic orbits (associated to fixed points) is 
studied via numeric integration of trimer motion equations thus revealing a diffused chaotic behavior 
(not excluding the presence of regular orbits), macroscopic effects of population-inversion and self- 
trapping. In particular, the behavior of orbits with initial conditions close to the dimerlike periodic 
orbits shows how the self-trapping effect of dimerlike integrable subregimes is destroyed by the 
presence of chaos. 

PACS numbers: 03.75.Fi, 05.45.-a, 03.65.Sq 



I. INTRODUCTION 

A remarkable progress in the experimental de- 
sign has been done since the first direct obser- 
vation of Bose-Einstein condensation in diluite 
atomic gas [[l]. One of most promising devel- 
opment concerns the construction of experimen- 
tal devices in which condensates, achieved within 
complex geometries, interact with each other giv- 
ing rise to quantum effects that are observable 
at the macroscopic level M -m. In this re- 
spect, one should recall, for example, the (super- 
fluid) boson Josephson-junction arrays obtained 
by means of optical lattices that trap weak inter- 
acting BECs in periodic arrays of potential wells 
||. In parallel with the experimental work, an 
increasing attention has been devoted theoret- 
ically to study the dynamical behavior of low- 
energy states in arrays of Bose-Einstein conden- 
sates (BEC) where the number of lattice sites 
(namely the potential wells occupied by the con- 
densate) is very large j(| 0. 

Opposite situations, corresponding to 'lattices' 
formed by two or three interacting wells, have 
been investigated as well in various recent pa- 
pers (see Refs. ^ -@)- In particular, the two- 
well system (dimer) has been analysed thorough/ 



from both the semiclassical (mean-field) and the 
urely quantum viewpoint in Refs. Jl3| , [l4 |] and 



||10| , [y| [T(|, respectively. Such investigations 
have revealed how the nontrivial structure of 
dimer phase space causes many significant phe- 
nomena such as the symmetry-breaking effect 
(issuing oscillation modes that are isoenergetic 
but inequivalent), the onset of 7r-phase oscil- 
lations, the self-trapping of boson populations, 
and, quantum-mechanically, the occurrence of 
(parameter-dependent, nondegenerate) doublets 
in the energy spectrum entailing periodic self- 
trapping. 

The connection between the quantum and the 
semiclassical picture of many interacting bosonic 
wells has been illustrated in Ref. M] and, 
in view of the closed link between an array of 
interacting BECs and the Bose-Hubbard (BH) 
model 0, in Ref. Q within the BH 

model theory. We wish to observe that the semi- 
classical approach (corresponding to describing 
condensates within the Bogoliubov approxima- 
tion) is appropriate for interacting wells with 
macroscopic boson populations. However, the 
systems recently obtained, where condensate is 
distributed among many wells are able to pro- 
vide mesoscopic numbers of bosons per wells Q] , 
might be best modeled by using the space-mode 
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approximation ||, [17| stemmed from the second- 
quantized boson field theory. 

Compared with the dimer nonlinear behavior, 
indeed the system of three interacting bosonic 
wells (trimer model) exhibits a new phenomenol- 
ogy and provides a novel research topic, both 
theoretically and experimentally. In fact, despite 
the simplicity of the model, the trimer dynam- 
ics is affected by a strong inner instability lead- 
ing to the chaos onset. This feature originates 
from the combination of the model nonlinearity 
with the nonintegrable nature that distinguishes 
trimer dynamics within extended regions of its 
phase space. In respect of this, the solutions of 
the trimer dynamical equations [displaying non- 
linear periodic oscillations (dimerlike orbits) with 
possible self-trapping effect] recognized in Refs. 
(O , fl2| must be viewed as special subregimes in 
which dynamics is integrable. Concerning the 
instability of trimer dynamics, we wish to re- 
call that the presence of homoclinic chaos has 
been revealed for the asymmetric trimer in Ref. 
pp[ |, while, at the quantum level, the survival of 
breather configurations has been investigated in 
Ref. ||. 

In the present paper we perform a systematic 
analysis of the trimer dynamics directed to as- 
certain that the nonintegrability generates chaos 
and that this dominates the interactions of three 
BECs modeled in a standard semiclassical pic- 
ture. A rich scenario of dynamical behaviors 
emerges from our analysis, which confirms the 
extremely structured character of trimer dynam- 
ics and represents the natural prosecution of the 
work of Ref. jl^] in which we focused our atten- 
tion on the special dimerlike (integrable) regime 
of trimer and on the self-trapping effect. 

Indeed the analysis we perform appears to be 
topical in relation to the study of the dynamics of 
solitons [Q and of vortices |2^] on one-dimesional 
chains of BECs as well as of experimental archi- 
tectures obtained recently The observation 
of such chain excitations, in fact, must take into 
account the possible destructive action of inner 
instabilities whose influence is clearly manifested 
in the trimer chain. We emphasize the fact that 
the trimer chain is the simplest possible situa- 
tion in which interacting BECs turns out to be 
governed by nonintegrable equations. The dimer 
dynamics, in fact, is completely integrable. 

The paper layout is the following. In Sec. II we 
review the derivation of the space-mode Hamil- 
tonian for three coupled wells from the quan- 



tum field theory of bosonic fluids in the dilute- 
gas approximation, and present the semiclassical 
picture that describes coupled boson wells with 
macroscopic populations. In Sec. Ill we identify 
the set of fixed points of the trimer Hamiltonian 
equations and show how such points are associ- 
ated to periodic solutions (representing collective 
modes) owing to the dynamically conserved total 
boson number N. Sec. IV is devoted to study 
the second variation (with N = const) of the 
energy function for such extremal configurations 
in order to make explicit their nature. In Sec. 
V we perform a dimensional reduction of trimer 
dynamics by defining a new set of canonical vari- 
ables whereby the constant of motion N is incor- 
porated in the dynamical equations. This paves 
the way to the implementation of the Poincare 
sections' method. The chaos onset in the trimer 
dynamics is investigated in Sec. VI. First, in a 
qualitative way, by constructing the Poincare sec- 
tions for trajectories whose initial conditions are 
chosen in proximity of the fixed-point configura- 
tions identified in Sec. III. Then, quantitatively, 
by measuring the maximum Lyapunov exponent 
of such trajectories. Sec. VII contains conclud- 
ing remarks and comments on future work. 



II. TRIMER DYNAMICS 

The model for a chain (or more complex struc- 
tures) of M interacting bosonic wells can be de- 
rived from the quantum field theory for boson 
fluids (with a nonlocal ip A term) by implement- 
ing the M-(space)mode approximation [Q. If 
the boson fluid is a diluite gas of TV interacting 
bosons trapped in an external potential V e then 
its dynamics is generated by the local boson-field 
Hamiltonian 0, ||, |§] 



H 



dvip (r) 
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where m is the boson mass, Uq = Airf^a/m takes 
into account the two-body interaction strength, 
and a is the s— wave scattering length. The 
field ip(r) [?/> + (r)] is the Heisenberg field oper- 
ator that destructs [creates] bosons at position 
r. The nonlinear term has been written in the 
usual normal ordered form. In order to work 
out the trimer Hamiltonian we state some as- 
sumptions: First, V e is a three-well symmetric- 
shaped potential. Second, the lowest energy level 
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H 3 {Z,Z*) := (Z\H 3 \Z) 



r 3 

^3 = 1 



U\z 3 \ 4 - v\ Zj \ 2 - ^ (z*z J+1 +c.c.) , (1) 



of each well (within the approximation of V e in man |TC 
terms of a single-well parabolic potential) must 
be well separated from the higher energy lev- 
els [jl7|. Third, the binary particle interactions 
is not strong enough to significatively change the 
latter assumptions. 

Some furher approximations are necessary to 
make explicit the ground-state structure in H. 
Let rj (i = 1, 2, 3) be the locations of the minima 
of V e and let Vj = V(r — rj) be the parabolic 
approximation to the potential in the jth min- 
imum, so that V e (r) ~ V(r — rj) when r ~ rj. 
Also, let us introduce the eigenstates Uj that rep- 
resent the normalized single-boson ground-states 
with energy Eq of the local parabolic potential 

Vj. These states are only approximatively or- Such equations for Zj (notice that Zj = (Z\a,j\Z) 



(j = 1 = 4 on the trimer chain) with the equa- 
tions 

ihzi = (2U\zi\ 2 - v)zi - y(z 2 + z 3 ) , 
ifiz 2 = (2U\z 2 \ 2 - v )z 2 -^{z 3 + z 1 ) , (2) 
ihz 3 = (2U\z 3 \ 2 - v)z 3 - y(zi + z 2 ) ■ 



thogonal because of J d 3 rv,jUk — 5jk + Rjk but 
the residue Rjk is a quantity exponentially sup- 
pressed depending on the overlap between Uj and 
Uk- The analysis is restricted to those potentials 
for which Rjk << 1 thus making negligible such 
contributions. 

The picture of the system thus resulting sug- 
gests that the boson field operator can be ex- 
panded as 4>(r,t) = ^iUi(r)&i(t), where a,i(t) 
[af(t)] is the annihilation [creation] boson oper- 
ator (associated to the space- mode state Uj) that 
satisfy the commutation relation [dj(i), %(t)] = 
5ik ■ By substituting this expression in the many- 
body hamiltonian one can obtain, to the low- 
est order in the overlap between the single-well 



modes, the quantum trimer Hamiltonian pO , 11 



H 3 = Yl i \U{n i -l)-v]n i --^T,(i t j) (a^o,- +a ] a l j , 

where the site index i,j = 1,2, 3, v = —Eq, and 
the operators rii = a\a,i count the number of 
bosons at site i. Also, in H 3 



T = 2 d 3 ruj[Vj - V e ]u j±1 



U = ^ /rf 3 r| Sl 



represent the (interwell) hopping amplitude and 
the strength of the Coulomb on-site repulsion, 
respectively [|[ p7[ . 

The quantum dynamics involved by Hamilto- 
nian H 3 can be cast in a classical form by repre- 
senting the system quantum state through a trial 
state \Z) written in terms of Glauber's states \zi) 
(defined by a,i\zi) = Zi\zi)). By implementing 
the time-dependent variational principle (TDVP) 
and the procedure discussed in Refs. Jn|, |2j| on 
\Z) = Hi\zi), one obtains the effective Hamilto- 



and z* = (Z\a^\Z)) can be calculated from H 3 
via the Poisson brackets {z^.,Zj} = iSk,j/h fur- 
nished by the TDVP method. Those for z* are 
easily obtained by complex conjugation. Various 
aspects concerning the special dimeric subregime 
of Eqs. (^|), in which trimer dynamics is inte- 
grable, have been studied in Refs. [O, Jl2j|. 



III. FIXED POINTS AND PERIODIC 
ORBITS 

The distinctive features that characterize the 
dynamics of a given Hamiltonian system can be 
recognized by exploring the structure of its phase 
space. The first step to do this is to locate the 
fixed points. The latter are derived in the present 
section based on Eqs. (||), whereas the nature of 
such points is studied in the next section. 

The fixed-point equations for the trimer are 

. Since the 



obtained by setting Zj = in Eqs. 
trimer dynamics is constrained by the constant of 
motion N — |zi| 2 + |z2| 2 + |2:3| 2 , one must incorpo- 
rate explicitly the restriction to the phase-space 
submanifold defined by N = const through a La- 
grange multiplier \. This just requires that one 
considers the variations of Ti 3 — x N in place of 
that of H. 3 . The resulting equations are 



(3) 



0=l2U\ Zj \ 2 -n+^ Zj -^Z 



where j ~ 1, 2, 3, and Z := z± + z 2 + z 3 , fi := x+w 
have been introduced . Therefore, any fixed- 
point, namely any vector (zi) = (rjx, r] 2 , i] 3 ) that 
satisfies Eqs. (||), provides, at the same time, 
a dynamically active solution of Eqs. (||) rep- 
resented by Zj(t) = r/j exp[ixt/K\. Despite its 
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time dependence, such a periodic orbit represents 
a (one-dimensional) extremal configuration of Tl 3 
on the hypersurface N — const. The solutions of 
Eqs. ^ can be grouped in various classes the 
first of which is represented by the case 



Z = Ej Zj = . 



The rem aining cases are obtained from Eqs. 
( A3)-( A4), where Zi can be replaced with the real 
quantities Xi (see appendix [A|). Such cases are 
given by 



x\ = x 2 = x 3 ^ , 
xi = —x 2 ^ x 3 = , 
x\ = x 2 ^ x 3 ^ . 
Such configurations are discussed below 

A. Ground-state configurations 



When xi = X2 = X3, Eqs . ( |A3| ) and (A4) are sat- 
isfied. Based on Eq. ( A2) combined with the con- 
served quantity N, one finds that Xj — ±y / N/3 
which provides 



z 4 = v/iV/3 exp(i$) , (4 = 1,2,3), (4) 

and 3/i = 2U N—3T. The energy of such configu- 
rations (they are shown to represent the ground- 
state in Sec. Ill) is given by 



E gs = \UN 2 - TN. 



(5) 



The ground-state phase $ is arbitrary since it 
represents a symmetry of the model. Also, the 
fact that Zj — \zj\ exp(i<f>j) have the same phase 
4>j = <f> reproduces the symmetry breaking phe- 
nomenon that distinguishes the minimum energy 
state and, particularly, the vanishing of the phase 
difference (between closed points) in superfluid 
media. Such two features, which are known to 
characterize the superfluid ground-state of BH 
lattice model, naturally extend to the trimer 
model of condensates in our semiclassical pic- 
ture. As regards Eqs. (||), they are solved by 
Zj(t) = y/N/3 exp[z($ + xt/K)], with x = l*~ «. 
that describes the ground-state collective mode. 



B. Vortexlike configurations 

The situation in which Z = leads to a special 
configuration. In this case Eqs. (0) reduce to = 



[2U\zi 



T/2) Zj implying, in turn, that 



2ju-T 
4U 



4NU+3T 



where the value of u is derived from Si|zi| 2 = N. 
As a consequence of the independence of \zj\ on 
the site index, the condition Z = is realized 
only if the phases of Zj = \zj\exp(i9j) are such 
that 9j(k) = 2njk/3+$ where k = 1,2 and $ is 
an arbitrary phase. This configuration represents 
a particular case of the vortex state discussed 
for the Bose-Hubbard model on a M-well chain 
lattice (M > 2) in Refs. || and |To|. The energy 
associated to the vortex states 



is given by 

E v 



UN 2 
3 



TNcos(^k) = 



TN 
2 



(6) 



(7) 



while dynamics issued from Eq. (g) is described 
by the solution Zj{t) = Zj(k) exp[ii(/i — v)/%\. 



Configurations with a single depleted 
well 



These configurations are characterized by the 
presence of a single depleted well (SDW). With- 
out losing generality, one can choose the second 
well, so that one has x±,X3 7^ 0, a nd x i = 0. This 
case can be faced based on Eqs. (Al): to satisfy 
E%t(X) = one must impose x 3 — —Xi, which 
entails, in turn, that the equation -Ei(x) = is 
equivalent to the equation -E^x) = 0. The lat- 



ter, together with the constraint TV 



<<3> 



implies that x\ = ±y/N/2 
provides the fixed points 



-X3,X2 = 0, which 



= ^Nj2c^, z 2 =0,z 3 = ^Wj2e^+*\ 



with fx = 2Ux\ + T/2. Their energy is 



E dw = \UN 2 + \TN. 



(8) 



(9) 



Permutations of site indices j of Zj furnish other 
five fixed points of the same type. In view of state 
(pi), it is worth noting that SDW configurations 
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have the same structure of the tt- states occurring 
in the dimer dynamics |L3[ , where the phases of 
each well keep a constant phase difference 7r in 
the course of time evolution. Solutions Zi(t) = 
ZjiO) exp[it(fi — v)/h] -with Zi(0) given by Eqs. 
((|)- exhibit this property. 



D. Dimerlike configurations 

For such states the variables Xj (recall that x 2 
is the boson number of the j-th well) satisfy the 
condition X4 = Xj =/= Xk- Three cases are thus 
obtained through index permutation. In general, 
since three fixed points turn out to be associ- 
ated with each dimerlike configuration (namely 
to each choice Xi = Xj ^ Xk), nine fixed points 
are finally found in the dimcric class. 

To deal with an explicit case, we shall con- 
sider the fixed points of the case x\ = x 2 ^ x 3 . 
Owing to its complexity, their derivation is de- 
scribed in appendix |B| Here we discuss the re- 
sults obtained. Fixed points depend on the pa- 
rameter t := T/NU and are representable as 
points (xi,X2,xa) G R 3 on a sphere due to the 
constraint N — Sjxf . Their expressions read 

Ai := (0,1,0,1, -(ai/\ai\)y/N - 2a 2 ) , 

A 2 := (a 2 , o 2 , -(a 2 /\a 2 \)^N ~2a 2 ^j , (10) 



A, 



where 



a 3 = ± 



03, 03, 



Np 2 



-{a 3 /\a 3 \)^N -2a 2 



l + 2p 2 



± 



N 



(11) 



and v — 1,2. Parameters q v and p are defined as 
implict fun ctio ns of parameter r through systems 
(B3) and (B4), respectively. In particular, the 
cubic equation 



r(2 + q 2 )(2 + q)+Aq(l + q) =0, 



(12) 



derived from Eq. (B3), with q £ [—1,0] fur- 
nishes the t functions Oi(r), g 2 (r) correspond- 
ing to the real roots of Eq. (|l2|). Notice that 
— 1 < <7i( r ) < <Z2( r ) < for r < r* , whereas, 
for t > t* , there are no solutions. The value t* , 
where qi(r*) = q 2 (T*) (and Ai = A 2 ), is calcu- 
lated in appendix pi. Instead, system (B4), with 



p 6 [—1/2,0], always exhibits a single solution 
p(r) that is carried out from the equation 

r(l + 2p 2 )(l + 2p)+4p{l+p) = 0. (13) 

The dependence on physical parameters T, U, N 
of Ai (qi ) , A 2 {q 2 ), A3 (p) by means of parameter 
r is thus established. 

Points Ai found in this way generate, by vary- 
ing r, three curves on the sphere with N — const 
[in view of the sign ± in Eqs. ( [Tl| ) they actually 
are six]. These become eighteen when consider- 
ing the fixed points generated by index permu- 
tations. This process is described in appendix 
[b|, where the actual number of dimerlike fixed 
points is shown to reduce to twelve. Such curves 
[parametrized by r via p{r) and <7„(t)] can be 
proven to never intersect with one another ex- 
cept for the special case r = r* where A\ = A 2 . 
This coalescence effect is discussed below. 

If the values of q v and p for some given r < t* 
are carried out explicitly, the energy for Ai 



E d = U[N 2 +6at-4Na 2 }-T[a 2 -2a^ N - 2a?] , 

_ (14) 
(i = 1, 2, 3) is obtained via formulas (|ll|). 

We conclude illustrating the physical situa- 
tions that correspond to configurations Ai, A 2 
and A3 when r changes. Also, we compare them 
to the pure-dimer scenario fl(i| , fL3| . Let us start 
with t — * 0. One has qi = —1 and q 2 = p = 
that entail 

Ai(-1) :=±(VA73, , 

A 2 (0):=( V ^V72, V ^V72, 0), (15) 



A 3 (0) 



respectively. By increasing r, Ai(qi) and A 2 (q 2 ) 
get closer and closer (Ai — A 2 for r — > r*). 
When t > r* only the fixed point A 3 (p) survives. 
In particular, r — > 00 implies p = —1/2, so that 
[from the third equation of Eqs. (fiO[) ] one obtains 

A 3 (- 1 -) = ±(^N/6, y /N/6,-^2N/3y (16) 

Since A 3 (p) [as well as B 3 (p) , C 3 (pV, obtained via 
index permutation (see appandix |Bj)] is shown to 
be a maximum in Sec. IV, indeed A 3 (p) appears 
to be comparable with the r-dependent maxi- 
mum of the (pure) dimer model [p0[, where a 



0, 0, 
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unique well ends up to be filled when t — > 0. In- 
stead, when t is increased, no merging of A 3 (j>) 
with other maxima [e.g., B^{p), C^(p)] happens 
as that observed in the (pure) dimer model. In 
this model, in fact, a (macroscopic) coalescence 
effect takes place [sec Rcf. |to[ since two 
maxima merge in a unique one when r > 1 [the 
opposite effect (bifurcation) occurs for r < 1]. 
As shown in Appendix A, such effects involving 
maxima pairs do not distinguish trimer dynam- 
ics. 

A different macroscopic effect happens, how- 
ever, in the trimer phase space. This is caused by 
the merging of distinct saddle points (e. g., Ai, 
A2) discussed above, that disapper for r > r*. 
Since chaotic behavior develops around saddle 
points, their coalescence should entail an evident 
local regularization of dynamical behavior. 

As a final remark, we wish to observe that 
state Ai(— 1) exhibits the same per-well boson 
distribution of the ground-state. Such two states 
differ uniquely owing to the phase of the third 
well. Similarly, state Ai(0) and the 7r-like state 
(y/N/2,-y/N/2,0) just differ owing to the op- 
posite phase of the second well. Despite their 
identical boson distribution, such situations will 
exhibit very different behavior thus confirming 
the profound influence of the interwell phase dif- 
ferences in distinguishing dynamical states. 



IV. CHARACTER OF FIXED POINTS 

In this section we consider the stability char- 
acter of fixed points just identified. Such a char- 
acter is recognized by studying the second varia- 
tion of the energy function on the hyper-surface 
defined by N — const. Explicitly, one should 
analyze the signature of the quadratic form as- 
sociated to the Hessian of H.3 (with N = const) 
in each fixed point. Being this process rather 
technical, we mainly develop it in Appendix [TJ. 
Below, after showing how the separation of H.3 
in two independent (local) subhamiltonians sim- 
plifies remarkably the stability analysis, we sum- 
marize the results obtained. 

To simplify H.3, it is advantageous to introduce 



the local variables & 



Vj, where Vj are 



the coordinates of some given fixed point. Then, 
neglecting third and fourth order terms, H.3 takes 
the form 

H 3 = H 3 (v) + Sf =1 (2U\ V] \ 2 -n-T) l&l 2 



+ >:; ,/•;/■-, + v 3 z*) 2 + Js^|& - ai 2 , (17) 

which undergoes the further semplification = 
v* — > Xi when one recalls that the phase factor 
of Vi € C is a constant that can be absorbed 
by £j. This fact allows us to separate H.3 in a 
q-dependent part and a p-dependent part, with 
£j := Qj + ipj. By making explicit the latter 
definition in TL3, we find 

H 3 =H + h(q ; 6U, T, ft) + h(p ; 2U, T, ft) , (18) 

where Hq := Hsiy), and 

h( q ■ en, t, n) := Y^ M *)a liU , (19) 



h(p ; 2U, T, ft) := ^(M p ) y p iPj . (20) 
ij 

Dynamical matrices M q and M p are defined as 

(21) 



T 

'2 



Ai 1 1 
1 A 2 1 
1 1 A 3 



T 
2 



<5i 1 1 
1 S 2 1 
1 1 83 



(22) 



with Aj := 2(ft - 6Ux?)/T and Sj 



2( M 



2Uxj)/T. Hence the diagonalization of the (lo- 
cal) quadratic form associated to H.3 can be per- 
formed in a separate way on ( JT9| ) and (20). Fur- 
ther simplifications come from the fact that ma- 
trix M p is proven to always have a vanishing 
eigenvalue (see Appendix |cj), while, due to the 
conserved quantity TV = Ei|zi| 2 = const, the in- 
duced local constraint qiXi + q 2 x 2 + (73X3 = 
makes h(q ; 2U, T, ft) dependent only on two vari- 
ables qi (see Appendix ^) . 

The stability character for ground-states, vor- 
tex states, SDW states, and dimerlike states is 
studied explicitly in Appendix [d| The calcula- 
tion of the H3 second variation and, when this 
is necessary, of Hessian eigenvalues provide the 
following scenario. 

- States with x\ — x 2 — X3 are energy minima; 

- Vortex configurations are saddle points; 

- SDW configurations are saddle points; 
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- Dimeric states exhibit two saddle points and 
one maximum point for < t < r* . 

In the last case the merging of the two saddle 
points to form a regular point is enacted for r — > 
t* . Thus, for r* < r, a single maximum survives. 

V. MAPPING OF DYNAMICS ON THE 
REDUCED PHASE SPACE 

We develop a both qualitative and quantita- 
tive analysis of the chaotic behavior of the trimcr 
based on the Poincare section method. After per- 
forming a qualitative study of the periodic-orbit 
instability we effect a quantitative analysis by 
measuring the very chaos' indicator, namely the 
Lyapunov exponent of every single orbit choosen 
near a period one. To this end, we introduce a 
more suitable coordinate system involving a sym- 
plectic reduction of dynamics. 

Hamiltonian (|l|) depends on three complex de- 
grees of freedom and commutes with the boson 
number N. This fact as well as the structure of 
coupling term in H.3 permit us to reduce from 
six to four the number of (real) canonical coordi- 
nates. The numerical integration of Eqs. (||) thus 
furnishes a system picture consisting of a trajec- 
tory in a four-dimensional {Ad) reduced phase- 
space ? € C 3 . In V 1 a Poincare section (PS) is 
the figure made by the points where a trajectory 
cut a 2d reference plane. The new set of canonical 
coordinates used to construct the PS embodies 
explicitly the conserved quantity E| =1 |^j| 2 = N. 
Complex coordinates — ^/n~i exp(iOi) are re- 
placed by 

01 = 2 - 6>i , £1 = (n 2 + n 3 - m)/N 
4>2 = 03 - 2 , £2 = ("3 - ni - n 2 )/N 
V> = (0 3 + 0i)/2 , N = m + n 2 + n 3 , 

(23) 

which obey the canonical Poisson brackets 

= -2<y tf /jv , {&,jv} = o, 

(24) 

ty,&} = o, {^,n} = -i. 

With such new variables Q H.3 becomes 

£ = 7jff2- W (£l> 01,02) = 

£r; + £ 2 2 -£i6-£i+6 + £o- 



r v / (i-ei)(ei-6)(i + 6)x 

COS 01 COS 2 COS 012 

_V^i + &J v / (Wi7 V(£i-6)J 

in which 0i 2 := 0i + 02, £0 = 1 — ^ n /{UN) 
(the associated Hamiltonian equations are con- 
tained in Appendix [|]) . In terms of coordinates 
(£1, £2, 0i, 02) the ground-state configuration (Q) 
and vortexlike fixed points @ correspond to 

(|,-|, 0, 0), (l-|,§7rfc,§7rfc) , 

(fc = 1,2), respectively. Dimeric fixed points 
[consider, e. g., A% in Eq. (|i~fj|)] are given by 

(a,2a-l,0,7r), 

where £1 = 1 — 2cl\{t)/N and a\{r) is defined by 
Eqs. (pT|), whereas SDW states can be expressed 
as (0, 0, <p, 7T — y>) for n 2 = 0. 

Operationally, the motion equations (|^) are 
numerically integrated by using a first-order bi- 
lateral sympletic scheme; the algorithm precision 
is checked by monitoring the conserved quanti- 
ties, that is, the system energy and total number 
of bosons. Trajectories can be traced in the phase 
space V in terms of £1, £2, 01, and 02. For any 
given value of the reduced energy, Hamiltonian 
( p5| ) defines a 3d hypersurface in V. The 2d sur- 
face used to construct a PS then is obtained by 
firming the value of £2 to a constant. Hence, PS 
is made by coordinates (£1, 0i) of the points' set 
in which trajectories cut the selected 2d surface. 

A. Discussion of numerical results 

We present the results of the numerical analy- 
sis aimed at investigating the phenomenology of 
the trimer dynamics in proximity of fixed points 
-these identify with the periodic orbits that sta- 
tionarize H.3 with N = const- calculated in Sec. 
III. The Hamiltonian parameters and the total 
boson number chosen for the numerical simula- 
tions are 

U= T= 1, v = 0, N = 10, 

respectively. Simulations have been carried out 
by using an integration time step of order 1 x 
10~ 4 , while the total number of time steps em- 
ployed in constructing each orbit is of order 2 28 . 
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For each periodic orbit a PS has been selected 
by setting £ 2 = const and d^/dt > 0; on this 
section we have considered samples of about 100 
initial conditions (IC). Also, the maximum Lya- 
punov exponent (MLE) has been measured for 
each trajectorie of the samples associated to the 
extremal periodic orbits. 

Ground-state. The PS of this case is fixed by 
the condition £2 = —1/3. About 100 initial con- 
ditions have been chooscn in the £1 — <j>\ plane, 
placed close to the point (£i,0i) = (1/3,0). 
The energy of the corresponding trajectories is 
E ~ 0A7(UN 2 /2). Fig. | shows that all the tra- 
jectories are regular in the phase space V . Their 
regular character appears to be consistent with 
the periodic character shown by the oscillations 
of populations rii — \zi\ 2 (see Fig. ||) of each 
condensate. 

Vortexlike initial conditions. The choice £2 = 
-0.3, E ~ 0.77([/A 2 /2) distinguishes the PS of 
the present case. The initial conditions for the 
trajectories have been chosen close to the fixed 
points (£i,0i) = (l/3,2fc7r/3), k = 1,2 (vor- 
tex state). Since no difference distinguishes the 
PS with k = 1 and that with k = 2, we re- 
strict our attention to k = 1. Fig. || shows the 
presence of both regular and chaotic trajectories 
[Fig. U supplies two examples, one for each or- 
bit type]. It shows as well that the PS points 
related to chaotic trajectories are distributed in 
a region well separated from that occupied by 
points generating regular orbits. In particular, 
the PS displayed in Fig. || [together with other 
PS involving slightly different £2 (— — 1/3)] sug- 
gests that the vortex-state fixed point is basically 
surrounded by chaotic orbits. In Fig. |3| to reach 
the nearest regular orbits starting from £1 = 1/3, 
4>\ = — 2n/3, a finite variation of both £1 and <f>\ 
is necessary. 

When this change is carried out the trimer- 
population oscillations change in a significant 
way. The time evolution of condensate popula- 
tions ni (i — 1,2,3) related to the nonchaotic 
orbit of Fig. |] is illustrated in Fig. ^ and con- 
firms its regular character. By considering IC 
closer and closer to the vortex state position such 
a character is progressively lost. This is shown 
in Fig. |^ that plots populations function 
of time, for the chaotic trajectory of Fig. |j. 

The regular orbit of Fig. [| involves an evident 
self-trapping effect provided their IC are enough 
far from (£i,0i) = (1/3, — 2n/3). This is clearly 



manifested in Fig. || where ni(t)'s oscillate in 
such a way that ri2(t) < rii(s), 71.3(5), Vi, Vs: a 
stable gap, in fact, separates ri2 from ri\, 77,3. On 
the contrary, no stable gap is involved, in general, 
by chaotic orbits (see Fig. |6[) that develop large 
oscillations on the whole range of n,-. The fact 
that, in average, n; ~ 1/3 is the only feature 
inherited by the vortex state. 

SDW-like intial conditions. In this case the 
choice £ 2 = -0.002 and E ~ 1.09([/iV 2 /2) fixes 
the PS that is represented in Fig. [?]. The zoom 
of the section reveals regular trajectories in V 
placed near SDW fixed points the latter being 
characterized by £2 = = £1, 4>i = <p, (f>2 = tt— <f, 
where we have set (if can be choosen arbitrar- 
ily) ip — 0.57 7T. Fig. H describes the PS of a 
chaotic orbit chosen among those of Fig. [jj The 
PS points are distributed in two, well separated, 
basins in a quite evident way. The interpretation 
of such an effect is the following: after recalling 
that setting £2 — implies that U3 — 1/2, one 
deduces that, concerning the points of the PS, 
the values allowed for n\ are either n\ ~ 1/2 or 
tlx — 0, which involves either — or 77,2 — 1/2, 
respectively. 

One therefore recognizes the presence of an in- 
version population phenomenon beetwen n\ and 
?i2. Interestingly, no intermediate values seems 
to be permitted. The corresponding scenario is 
given, on a shorter time interval, in Fig. || where 
the nonperiodic oscillations of populations fti(t), 
n2(t), and nz(t) are compared and a population 
inversion effect involving Tlx, ni gets going. The 
population oscillations referred to a regular or- 
bit of those contained in the zoom of Fig. are 



shown in Fig. 10 



Dimerlike initial conditions. As proven in Sees. 
Ill and IV, the fixed points of this case con- 
sist of two saddle points and a maximum. The 
conditions £ 2 = -0.295, E ~ 0.73(UN 2 /2) and 
£2 = -0.005, E ~ 0.91(C/iV 2 /2) firm the PS as- 
sociated to the first saddle (see Fig. |ll|) with 
coordinates £1 = —£2 = 0.295, <j)\ = —02 = 7T, 
and to the second saddle (see Fig. |l^) with co- 
ordinates £1 = —£2 = 0.005, 4>i — —02 = 7T, 
respectively. In both the cases, the PS's exhibit 
both regular and chaotic trajectories. Concern- 
ing Fig. [n] (first saddle point), the coexistence of 
such regimes is confirmed in Fig. [l^, where the 
PS of a regular trajectory is compared with the 
PS of a chaotic one. Figs. [l6] shows analogous 
quantities referred to the second saddle point. It 
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is worth noting that in the latter case the neigh- 
bourhood of the saddle point is characterized IC 
issuing regular orbits, whereas the first saddle is 
surrounded by IC generating chaotic motions. 

In Fig. [l5|, regular orbits reside either on the 
bottom or on the top of the figure, but dynamics 
never connects the top orbits with the bottom 
ones. This feature is confirmed by the time be- 
havior (see Fig. [l7]) of rii(t) related to the regular 
trajectory of Fig. [l6| In Fig. [l?], the presence 
of the gap beetwen (— 0.5) and m, ^3 (~ 5) 
indicates a macroscopic self-trapping effect. The 
latter differs from the dimeric self-trapping re- 
viewed in Sec. II-A, where Z\ — 2:3 — * n\ = TI3, 
in that n\, 77.3 develop independent oscillations. 
The scenario just described no longer holds for 
the chaotic orbits: Fig. [l^ shows an intermittent 
effect of population inversion beetwen n\ and ri2. 
This reflects the fact the points of the chaotic or- 
bit of Fig. |l6| are distributed intermittently both 
in the higher and in the lower part of the PS. 

The trajectories near the maximum (we con- 
sider the case £1 = 0.999, £ 2 = 0.998, fa = 
—fa = 1")) appear to be regular, as one can 
deduce from Fig. [1^ that shows PS close to 
a maxium fixed point with £2 = 0.8 and E ~ 
1.7(UN 2 /2). In this case, populations n\{t), 
riiif), and n^{t) display in Fig. ^ a periodic 
effect of population-inversion involving n\, 712, 
whereas the fact that m, n-i << 713 entails an 
evident self-trapping phenomenon. 

The averages of the MLEs calculated for the 
chaotic orbits are 0.74 (Fig. 3), 0.42 (Fig. @), 
0.29 (Fig. [n|), and 0.96 (Fig. [jj). The evalua- 
tion of the MLE for the regular trajectories de- 
scribed in the previous examples (as well as for 
those of Fig. [I] and Fig. [l9|) exhibit the expected 
decreasing behavior thus suggesting that a weak 
stochasticity occurrs on such trajectories. 



VI. CONCLUSIONS 

In this paper we have focused our attention 
on the structure of the phase space of 'classic' 
trimer, that is the mean-field form of the model 
describing three interacting BECs. Our analysis 
puts in light, on the one hand, the remarkable 
complexity that characterizes the trimer dynam- 
ics (by comparison with the integrable dynamics 
of the dimer system), on the other hand, the phe- 
nomena that are expected to characterize the M- 
well chain of interacting condensates. In view of 



the recent experimental results, the phenomenol- 
ogy of this system seems more and more viable to 
experimental observations. Trimer dynamics has 
been investigated within a semiclassical Hamilto- 
nian picture, reviewed in Sec. II and developed in 
previous papers, based on a coherent-state rep- 
resentation of the trimer quantum state. 

The identification of the set of fixed points 
of trimer Hamiltonian equations and the fact 
that such points are associated to periodic so- 
lutions (collective modes) of several types rep- 
resents the initial, central result of our paper. 
The presence of the constraint N = const en- 
tails that the states that stationarize the Hamil- 
tonian are not isolated points but periodic orbits 
(one-dimensional manifolds). The solutions thus 
found enlarge the set of exact solutions Q, [l^] 
pertaining the dimerlike integrable subregimes 
of trimer and exhibiting a parameter-dependent 
self-trapping effect. 

Based on the second variation (with N — 
const) of the energy function around its fixed 
points, the character of the latter has been recog- 
nized in Sec. IV and appendix |5| thereby reveal- 
ing the presence of several saddle points and max- 
ima, in addition to the expected ground-state. 
Numerical simulations and the PS method have 
furnished a wide scenario of trimer dynamical be- 
haviors whose possible chaotic character has been 
tested by measuring the maximum Lyapunov ex- 
ponent. We summarize the results of our dynam- 
ical simulations. 

i) The orbits that have ICs close to ground-state 
exhibit a regular behavior with periodic oscilla- 
tions of populations rij. 

ii) On the contrary, orbits with vortexlike IC 
(namely based at points close to vortex fixed 
points) are, in general, chaotic. Regular orbits, 
however, are found at sufficient distance from 
vortex fixed points. For such orbits a stable gap 
separates the oscillations of rii from those of rij , 
n k (i) k i) thus generating self-trapping; rij, n^ 
show independent oscillations (in Sec. VI con- 
figurations with i — 2, j — 1, k — 3 has been 
considered). The gap disappears for chaotic or- 
bits. 

iii) Orbits with SDW-like ICs also display both 
regular and chaotic behaviors but their ICs are 
not separated spatially. The regular orbits we 
have considered keep memory the IC since one 
of the three wells remains almost empty (pure 
dimer) while the other two undergo regular os- 
cillations. Such states identify essentially with 
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it- like states and manifest a stable character. An 
example of chaos emergence, which starts with 
a macroscopic population inversion entailing the 
filling of the (initially) depleted well, has been 
detected by assuming various SDW-like ICs in- 
volving a chaotic orbits. 

iv) Regular orbits generated by dimerlike ICs 
(related to the second saddle point) exhibit pe- 
riodic oscillations of rij's with an evident self- 
trapping. Chaotic orbits, instead, give rise to os- 
cillations displaying self-trapping on short time 
intervals and intermittent population-inversion 
effects. ICs close to the maxima further gener- 
ate regular orbits with self-trapping. The corre- 
sponding states display the presence of a unique 
almost filled well 

The scenario just depicted supplies a rich ac- 
count of properties, behaviors and possible ob- 
servable effects issuing from trimer dynamics and 
suggests promising future developments. An ef- 
fect that might have a macroscopic character is 
related to the expected chaos suppression (see 
Sec. HID) caused by the A\-A?, coalescence for 
r = t*. The basic configurations (ground-state, 
vortex states, SDW states, and dimerlike states 
involving both saddle and maximum points) rec- 
ognized in the present work, and the complex- 
ity of dynamical regimes, both chaotic and reg- 
ular, that develop in their neighbourhoods de- 
serve further investigations in two directions at 
least. First, classically, one should carry out a 
systematic study (requiring huge computational 
resources) of long-time behavior of dynamical 
states of interest to disclose further macroscopic 
effects. For the same reason, a larger number of 
ICs (together with the trajectories thus issued) 
should be considered near fixed-point configura- 
tions. 

We point out the fact that predictions on 
the dynamics of phases fa should be important 
in relation to phase-interference experiments |S . 
This aspect, which has not been deepened in the 
present work, requires a separate analysis and 
further numerical study directed to detect phe- 
nomena exhibiting phase coherence and their sta- 
bility in proximity of states endowed with or- 
dered phase configurations such as vortex states 
(fa = 2-Kk/3, k = 1, 2), SDW states (fa - fa = tt, 
i ^ j, \zk\ =0 with i,j ^ k), and dimer configu- 
rations (fa = fa, Zk Zi = Zj with i,j k). 

Second, in view of the possibility of realizing 
systems with small per- well populations, the pure 
quantum approach to trimer dynamics (along the 



same lines of previous work directed to study the 
spectral properties of dimer) seems to be quite 
natural. The study of quantum trimer might 
put in evidence unexpected effects caused by the 
competition of chaotic (classical) behavior and 
integrable (quantum) behavior on the border- 
line of appropriate mesoscopic regimes where the 
transition from quantum to classical dynamics 
takes place. 
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APPENDIX A: DERIVATION OF FIXED 
POINTS 

Eqs. ^ can be simplified by noting their in- 
variance under the global symmetry transforma- 
tion Z£ — > Z£ exp(i<i>) and the fact that z^/Z £ R 
(whenever Z ^ 0) with t = 1,2,3. Then one can 
set Z£ = X£ exp(i$), where the xi are real num- 
bers and $ is an arbitrary phase, thus reducing 
Eqs. (||) to a system of the three real equations 

Ej (x) := [2Ux] - n + ~) Xj - ~X = , (Al) 

with j = 1, 2, 3 and X := x\ + X2 + x 3 . When 
the condition Xi ^ for i = 1, 2, 3 is imposed, 
the latter equations can be recast in terms of an 
equivalent system of three equations one of which 
fixes the Lagrange multiplier fi, while the other 
two, now formulated in a /i-independent form, 
determine x\, X2, X3 thanks to the further con- 
dition N — const. In fact, the sum of the quan- 
tities Ej (x) / Xj can be set equal to zero provided 
Xj 7^ thus giving the equation 

3/i = 2C/iV-|s i ^^. (A2) 
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Moreover, from E\{x)/x\ — E 2 (x)/x 2 — and 
E 3 (x) jx 3 — E 2 (x) jx 2 — one obtains 



= (x 2 - xi) 



= (x 2 - x 3 ) 



2U(x 2 +x 1 ) + 



2U(x 2 +x 3 ) 



TX 



2x\X 2 



TX 



(A3) 



(A4) 



2x 3 x 2 

completed by the condition N — x\ + x\ + x 2 . 



APPENDIX B: DERIVATION OF FIXED 
POINTS OF DIMER TYPE 

Because of the identification x± — xj__ (character- 
izing dimerlike fixed points), Eqs. ( [A3] ) and (A4) 
become a unique equation that can be written as 



2U{x\ 



Xn = — 



1 



2-3 _ 2 Xl 



X\ 



x?, 



(Bl) 



Such an expression suggests two possible ways to 
parametrize x\, x 3 



(B2) 



xi = y'NRcha , x 3 — v NR sh a 
xi=VNRshf3, x 3 = y/NRchp 



both allowing for the elimination of N in the con- 
straint 2x\ + x\ = N. They also provide two 
independent class of solutions that make explicit 
the t hree roots involved by the cubic character of 
Eq. ( Bl ) . In the first case one finds 



9 1-9 
R 2 = ^ 



■r 



R 2 - T 
R "4 



2 

q- - 

q. 



(B3) 



[the first formula comes from the constraint on 
the total number of bosons, the second one comes 



from Eq. (Bl)], whereas the second choice gives 

„2 



l + 2p 2 



R 



2p- 



1 

1 - - 
P. 



(B4) 



where r := T/U N, q = tha, p = th/3, and a, /3 € 
R. System (B3) reduces to the cubic equation 



r(2 + q 2 )(2 + q)+Aq(l + q) = 0, 



(B5) 



which, provided q G [— 1,+1] in order to ensure 
the condition R 2 > 0, supplies either two or none 
solutions, depending on the fact that r < r», 
r > 7~*. By solving the system one finds that 



the two roots q v (r), v = 1, 2 range in [—1, 0] and 
fulfil the conditions 

-1 < qi < q2 < , for < r < r* , 

with qi(r) = q 2 (r) for r = r*. The parameter 
7 % i s identified by imposing, in addition to Eq. 
( |B5| ), the requirement that the two ^-dependent 
functions of Eqs. (|b|) (that is, the two right- 
hand sides) are tangent at some point 



-6q 



(2 + 9 2 ) 



_ T 
= 4 



1 



(B6) 



Equations (|B5j) and (|B6|), solved numerically, 
supply the value r* ~ 0.29718. 

Such a structure does not characterize system 



(B4), which always exhibits a single solution for 
some appropriate value in the sector p G [—1/2, 0] 
obtained from the equation 

r(l + 2p 2 )(l + 2p)+4p(l+p) = 0. (B7) 

In view of the restriction q, p < 0, from defini- 



tions (B2) one deduces that the fixed-point coor- 
dinates are such that x±x 3 < 0. The three solu- 
tions qi(r), q 2 (r), p(r) just obtained correspond, 
within the space of coordinates {(21, x 2 , x 3 )} = 
R 3 , to three vectors expre; 

Ax := (ai,ai,-(ai/|ai|VW-2a?) 
A 2 := (a 2 ,a 2 ,-(a 2 /\a^ 



as 


'N - 


2a\ 


<N- 


2a\ 



(B8) 



A 3 
where 



a 3 = ± 



Np 2 



1 + 2p 2 



± 



N 



ql 



(B9) 

with q„, p solving Eqs. (B5), (B7). One can eas- 
ily check that the two r-dependent curves A\ and 
A 2 can be seen as two branches of a unique curve 
based at the common point Ai(qi) = A 2 (q 2 ) for 
r = t* , where they join smoothly. 

Dimerlike fixed points become nine when con- 
sidering the further set of points generated by 
index permutations 

Si := (oi, -(ai/larD^N -2al,a 1 ), 



B 2 := (o a , -{a 2 /\a 2 \WN -2a 2 2l a 2 ), (B10) 



B 3 := {a 3 ,-{a 3 /\a 3 \)^/N -2a 2 ,a 3 ) 7 
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d := (-(oi/|oi|ViV -2o?,oi,oi), 



+ 2 - = 



(CI) 



C 2 := (~(a 2 /\a 2 \)y/N-2a 2 2 ,a 2 ,a 2 ), (Bll) 

C 3 := (-(03/1031^-24,03,03), 

related to the subcases 11 = 3:3^ x 2 and X3 = 
X2 ^ ij. Such nine curves become actually six. 
Our previous observation on considering Ai and 
A 2 as a unique curve, in fact, readily extends the 
other curves B\, B 2 , and C\. C 2 . Owing to the 
double choice ± in formula ( |ll| ) such curves are 
twelve. 

In order to visualize the dimerlike fixed points, 
one can plot their position vectors Ai,Bi,Ce 
t — 1,2,3 within the three-dimensional space by 
varying their parameters ag (£ — 1,2,3) in the 
appropriate range. As noted above, in the case 
r < t* , for q ranging in the intervall [—1,0], Eq. 



(B5) exh ibits two solutions, whereas, for any r, 
Eq. ( |B7| ) admits one solution with p 6 [—1/2, ] . 
The corresponding ranges of variation for a v 
{y = 1, 2) and 03 are 



N 1 



< 



N , , N 

j, o<H<V ¥ 



respectively. The representation of the position 
of dimeric fixed points on the sphere S^xf = N 
provides arcs that never intersect the one with 
the other when t ranges in [0, 00]. 

The double sign ± in parametrizations (p9) 
entails that each curve Pj (P = A,B, Cjis 
formed by two disjoint curves. For P 3 (P = 
A,B,C) one has 



where 5iS 2 d 3 + 2 — T,jSj, u pon setting Sj = — (x£ + 
Xk)/xj owing to Eqs. (Al), can be shown to 
vanish in virtue of the identity 



„„. X/ + Xk 

n* - — - = 2 + s* 



X£ + x k 



The superscript symbol * recalls that the indices 
£, k, and j must differ the one from the other. 
Hence, as a general result, the diagonalization 
of M p entails the presence of a zero eigenvalue 
consistent with the analysis of dynamics in the 
reduced phase space we develop in the sequel. In 
view of the invariance of the matrix trace, one 
also finds 



v 3 \ ■ — T v 3 
^i=i A J - 2 



Ai + A2 



T „ X t 

2 ^-j — ; 



Xk 



(= -3/i + 2UN) 



while the two roots 



K = 



-A± V(A) 2 -4(A 123 -3) (C2) 



where, A = S, Si and A123 := Sj 616263/6^, and 
v = 1, 2, [ Ai (A2) is associated with — (+) in the 
formula] , are derived from the quadratic equation 
that emerges from Eq. ( pl[ ) when removing a 
factor A. 
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a 3 e [-VN/6,0[, a 3 e]0,y/N/6], 

while for P x and P 2 (P = A, B, C) the two dis- 
joint branches are originated by the mappings 
a„ — * P v with 



[-y/N/2,-y/N/3] , a, € [VW, VWA 



APPENDIX C 

The diagonalization process of M p allows one to 
prove that one of the eigenvalues is always zero. 
^From the standard condition det(Af p — A) = 
one obtains the eigenvalue equation (upon intro- 
ducing A := 2A/T) 



This appendix is devoted to recognize, case by 
case, the minima, the maxima and the hyperbolic 
points within the four (class of) states identified 
as fixed points. 

Concerning the eigenvalues of matrix M q one 
must take into account the restriction on the dis- 
placements £j from V{ induced by the constraint 
N = S^Zjl 2 = const. After recalling that the 
phase of Vi can be absorbed, for each i, by £j 
due to its arbitrariness, the substitution Vi — > xi 
implies that 



N = |& 



Si [\Zi\ 2 +2qiXi 



which, in turn, entails J^i\^i\ 2 + 2qiXi 
-to first order- the plane equation 



0, namely 



A 3 + (E^-)A 2 + (M2 + S 2 S 3 + S 3 S 1 - 3)A+ 



q x Xx + q 2 x 2 + q 3 x 3 = 
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It represents the restriction on the displacements 
that variables £j's are allowed to effect. Substi- 
tuting qi with qi = — {x r q r + x s q a )/xi, where 
r, s 7^ i (and the choice of i depends in general 
on the condition Xi =/= 0) finally gives 



h(q ; 6f7, T, fj,) 



■-(6Uxl-n) + 



T(x r 



q r q s 



12Ux 2 -4 (a:? 



(Dl) 



Ground-state case 



These fixed points are characterized by the fact 
that Xj = ±y/N/3, for j = 1,2,3, and /i = 
2UN/3 — T. By inserting this solutions in Eq. 
(18) one obtains the Hamiltonian written as 



H 3 



E„ 



3 

E 

3=1 



r 4 



T)q 2 + Tpj 



1 X ^ 

~ 9" E + ' 



where E gs is the ground state energy defined pre- 
viously. Then, by taking the constraint ./V = 
^f=il z j| 2 m t° account, one obtains S| =1 [g| + 
p 2 ± 2yjN/3qj] = 0. For little displacements 
from qj = = pj the latter equation reduces 
to Tij =1 qj = 0, which implies that 



H, 



E n 



(-UN + 3T)(q 2 1 +q 2 2 +q 1 q 2 ) + 



7 2_, 



(D2) 



The eigenvalues of the Hessian corresponding to 
the g-dependent and the p-dependent part of 
Hz are {(3t + 8/3)UN/2, (9r + 8)UN/2} and 
{0,3T/2,3T/2}, respectively. They are positive, 
coherent with the fact this is a minimum. 



Vortex case 



The conditions \x l \ 2 = N/3 and (j, = (4NU + 
3T)/6 characterize vortex configurations. One 



thus finds 

h(p; 2U, T, n) = -T(pi +p 2 + P3) 2 /2 < , 



whereas from Eq. (Dl) one gets 
h(q;6U,T,fi) = WN{q\ 



4 + 9i©)/3 > 



whose eigenvalues are always postive. Vortex 
configurations are therefore saddle points. 



3. SDW case 

The previous analysis shows that the (fixed 
point) configurations in which one of the three 
well is depleted (e. g., well i = 2) is such that 
x 2 = 0, xi = — x 3 — ±yjN/2 (7r-state structure), 
and fi = NU + T/2. Site index permutations al- 
lows one to obtain two further, similar cases. In 
these points, Hamiltionan (J18|) can be written as 



Hz — Edv 



UN 



[2{q\ + qi) - q 2 ] -; 

[(E^) 2 + CE>) 5 



In this case, the constraint on the total number of 
particles supplies the constraint Y^=i(Pj +3j ) ^ 
V2N(qi — q 3 ) = 0, which reduces to q\ — q 3 ~ 
when qj,Pj — 0. Then H3 takes the form 



H 3 ~ E dw + UN(4q 2 - q 2 ) - UNp 2 



(2qi + q 2 f + (pi + p 2 + psf 



(D3) 



whose Hessian matrix is endowed with the eigen- 
values {(6-5t±^5(20 - 12r + 5r 2 ))C/V/4} and 
{0, -(2 + 3t ± V4-4t + 9t 2 )[/V/4} for the q- 
dependent part and the p-dependent part, re- 
spectively. The analysis of the signature of such 
eigenvalues leads to identify the fixed points of 
the empty- well case with saddle points. 



4. Dimeric case 

In the dimeric case the conditions on the co- 
ordinates are x :— xi — x 2 7^ £3 =: y (none 
vanishing). Furthermore, one has to impose the 
condition (2Ux 2 - fj,)xj = T/2 J2kjtj 

Xk on the 

Lagrange multiplier [i that becomes 



fi= (2UN -T(l + y/x + x/y))/3 



(D4) 
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Let us start by analyzing the Hessian eigenval- 
ues of Hamiltonian p sector namely of E q. (p0[). 
These eigenvalues are given by Eq. (C2) and one 
has that whenever the term a := (6162 + ^2^3 + 
S3S1 — 3) is greater then zero, these eingevalues 
are negative. To verify the condition a > one 
can proceed in the following way. First, one sub- 
stitutes in a the Sj written in term of thier defi- 
nitions; namely in term of [a, T, U and the fixed 
points coordinates x, y, Second, one writes u in 
terms of x, y and etc. [Eq. (p4|)] . Finally, one 
eliminates the dependence on x, y and T in a by 
choosing either the parametrization 



y/x = q, x = ±^N/{2 + q 2 ) 

combined with Eq. ( p35| ) to express the in- 
dependence of r, or the parametrization 



x/y=p, x = ±y/Np 2 /{l + 2p 2 ) 



combined with equation (|B7[ ) to express the p- 
dependence of r. The expression for a achieved 
in such a way depends on q or p, respectively. 
One can show that both the expressions are al- 
ways positive in the range of definition of q 
([—1,0]) and p ([—1/2,0]), which means that 
eigenvalues ( |C2| ) are all negative. 

As usual, for working out the Hessian eigen- 
values of Hamiltonian (091) related to the q part 
of the original one, it is necessary to take into 
account the constraint £i|zi| 2 = N = const. 
The latter, in the present case, becomes (73 = 
~(Qi + Q2)x/y. By means of this condition, one 
can reduce the dimension of the eigenvalues prob- 
lem related to Hamiltonian (|l^) from 9 to 4. The 
Hamiltonian h r {<li,(l2,U, jJb,T) thereby obtained 
can be further simplified through the substitu- 
tions fi -> u = [2UN - T(l + y/x + 2x/y)]/3 
and y/x = q or x/y = p, depending on the 
parametrization one adopts. With the first choice 
(y/x = q), and relying on Eq. (B5), one finds two 
Hessian eigevalues one of which is always posi- 
tive, whereas the other has an illdefined sign in 
the domain q 6 [—1,0]. By using the second 
parametrization x/y = p and Eq. (B7), both the 
eigenvalues thus obtained can be proven to be 
negative for p G [—1/2,0]. In summary, in the 
dimeric case, for < r < r*, one has two saddle 
points and one maximum point; for r* < r, in- 
stead, fixed points reduce to a single maximum 
point. 

APPENDIX E 

The fixed-point configurations corresponding to 
the change of Hamiltonian variables Zi, z* — ► £ a 



<f) a , N, tf> are obtained from the equations of mo- 
tion rewritten in terms of the new variables (see 
below). Coordinates transformation ( |23| ) can ex- 
hibit (isolated) singular points in which they are 
not invertible. Dimer configurations in the pres- 
ence of an empty well provides an explicit exam- 
ple where transformations ( |23| ) are illdefined. In 
fact, fixed points 



fa) = (V^V/2exp^,0, VA72exp^ + 7r)) 

correspond to the set {(0, 0,x,tt — x)l^x} ; i n the 
new description. As regards dynamical applica- 
tions, fortunately, this problem is bypassed be- 
cause the trajectories chosen in proximity of the 
periodic orbits associated to this kind of fixed 
points (as well as the PS used to study the dy- 
namics near the same fixed points) do not con- 
tain, by construction, such pathological points. 
Upon setting 0i 2 := </>i + 2 , s :— UNt, the 
equations of motion on the reduced phase space 
are given by 



-r = 1 - 2 6 

ds 



I 1 + 6 

6-6 



(1-2^+6) cos 0i /l + e 



S2 



V(i-6)(6-6) V 1 - ^ 



COS 012 



(El) 



U-02 * T 



' 1-6 

6-6 



COS 01 + 



(1-6 + 26) cos 4> 2 1-6 



V(6 -6)(i + 6) V 1 + £ 



COS 012 



(E2) 



ds 



1 = tV(1-6)(6-6) sin 0i + 



V(l -6X1 + 6) sin0 
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(E3) 



J=V(6 -6X1 + 6) sin 02 + 



V(l -6)(l + 6) sin 0i2 . (E4) 
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FIG. 1: Poincare section at E ~ 0A7(UN 2 /2) and 
£2 = —1/3. In the neighbourhood of the ground-state 
the phase-space trajectories are regular. 



FIG. 2: Tipical time evolution of the condensates' 
populations related to a motion with initial condi- 
tions close to a ground-state configuration. Solid 
line refers to ni(i), dashed line corresponds to ri2(£), 
713 (t) is represented by the dotted line. The dynam- 
ics appears to be periodic (within the simulation time 
scale) . 



FIG. 3: Poincare section at E ~ 0.77(UN 2 /2) and 
£2 = —0.3 for orbits close to a vortexlike fixed point. 
Even if some regular orbits are present, phase-space 
trajectories show a dominating chaotic character. 
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FIG. 4: Representation of a regular orbit (identified 
by the almost continuous line) and of a chaotic one 
chosen among those of the PS of Fig. |§| 

FIG. 5: Time evolution of condensate populations 
ni(£) (solid line), 712(f) (dashed line) and 713(f) (dot- 
ted line) for the regular trajectory of Fig. ^| 



FIG. 13: Temporal behavior of the condensates pop- 
ulations related to the regular trajectory of Fig. [12. 
The solid line, the dashed line, and the dotted line 
refer to 711(f), 712(f) an d 713(f), respectively. 



FIG. 6: Populations rii (solid line), 712 (dashed line) 
and 773 (dotted line) are plotted as a function of the 
time for the chaotic trajectory showed in Fig. W. 



FIG. 14: Condensates populations, as a function of 
the time, related to the chaotic trajectory of Fig. [l^. 
The solid line, the dashed line, and the dotted line 
refer to ni(f), 712(f) and 773(f), respectively. 



FIG. 7: Poincare section at E ~ 1.09{UN 2 /2) and 
£2 = —0.002 relative to orbits close to a SDW fixed 
point 



FIG. 

tories contained in the Poincare section of Fig. 



raje 
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FIG. 15: Poincare section firmed by the conditions 
£2 = -0.005 and E ~ 0.91(IW 2 /2), close to a saddle 
dimeric fixed point. 



FIG. 9: The figure plots the populations ni(f) (solid 
line), 712(f) (dashed line) and 773(f) (dotted line), as 
a function of the time, for the chaotic trajectory 
showed in Fig. H. 



FIG. 10: Populations of the three condensates: 7ii 
(solid line), 712 (dashed line) and 713 (dotted line) as 
a function of the time for a regular trajectory of the 
zoom showed in Fig. 0. 
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Poincare section firmed by the conditions 

2 , 



FIG 

£2 = -0.295 and E ~ 0.73(UN 2 /2), close to a 
dimeric (saddle) fixed point. 



FIG. 16: This figure shows a regular trajectory and 
a chaotic one chosen among those of Fig. hq. 



FIG. 17: Temporal behavior of the condensates pop- 
ulations related to the regular trajectory of Fig. hq. 
The solid line, the dashed line, and the dotted line 
trace Tii(f), 712(f) and 773(f), respectively. 



FIG. 12: This figure shows a regular trajectory and 
a chaotic one chosen among those of Fig. [ll]. 
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FIG. 18: Condensates populations, as a function of 
the time, related to the chaotic trajectory of Fig. |l6| 
The solid line, the dashed line, and the dotted line 
trace ni(t), U2(t) and ns(t) respectively. 



FIG. 19: Poincare section near a maximum at E ~ 
1.7(t/7V 2 /2) and £ 2 ^ 0.8. In the neighbourhood of 
the maxima the phase space trajectories are regular. 



FIG. 20: Time evolution of the condensates popula- 
tions related to a motion with initial conditions close 
to a maximum state configuration. Figure shows 
ni(t) (solid line), n2(t) (dashed line), and n^(t) (dot- 
ted line). The motion appears to be regular. 
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